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1.  INTRODUCTION 

This  report  documents  work  performed  on  AFOSR  Grant  90-0271  —  “An  Addidve 
Turbulent  Decomposition  of  the  Navier-Stokes  Equations  Implemented  on  Highly  Parallel 
Computer  Systems.”  The  work  was  performed  at  the  University  of  Kentucky  (UK),  and  (on  a 
subcontract  to  the  grant)  at  UCLA,  during  the  period  1  June  1990  to  31  May  1991.  A  repon 
on  the  first  four  months  or  effort  was  submitted  in  the  form  of  an  Annual  Repon  at  the  end  of 
October,  1990.  We  shall  begin  this  report  by  presenting  some  background  information, 
following  this  with  a  brief  summary  of  the  first  four  months  of  work  and  then  the  main  body 
of  the  report.  The  main  items  treated  in  this  report  will  represent  work  completed  during  the 
period  1  October  1990  to  31  May  1991;  they  include  our  continuing  work  on  the  transition  to 
turbulence  in  channel  flow,  the  beginnings  of  a  2-D  implementation  of  the  additive  turbulent 
decomposition  (ATD)  in  generalized  coordinates,  work  on  the  interaction  between  chemical 
kinetics  and  turbulence,  and  further  analyses  of  ATD  and  domain  decomposition.  We  shall 
end  with  a  summary  of  our  work  on  this  project  as  a  whole,  including  all  work  supported  by 
AFOSR  grants  90-0271  and  89-0281,  and  will  list  the  publications  and  presentations  which 
have  resulted  from  this  work 

1.1  Background 

Additive  turbulent  decomposition  (ATD)  is  an  approach  to  turbulence  simulation,  first 
proposed  by  McDonough,  et  al.  [1],  who  were  motivated  by  the  observation  that  none  of  the 
present  methods  for  turbulence  simulation  are  really  acceptable  for  the  broad  range  of 
turbulent  flows  that  must  be  analyzed  for  practical  engineering  problems.  The  classical 
methods  —  mixing  length,  k-e,  and  Reynolds  stress  models  —  are  efficient  enough,  but  are 
not  predictive.  Moreover,  if  the  Reynolds  averaging  procedure  on  which  these  methods  are 
based  corresponds  to  time  averaging,  die  results  will  not  even  be  consistent  with  the  Navier- 
Stokes  equations.  On  the  other  hand,  higher-level  techniques  such  as  large  eddy  simulation 
(LES)  and  direct  numerical  simulation  (DNS)  provide  results  that  are  generally  consistent 
with  solutions  to  the  Navier-Stokes  equations,  and  are  predictive  or  nearly  so,  but  they  are 
not  efficient  enough  for  application  to  high  Reynolds- number  (Re)  engineering  problems  on 
present  or  forseeable  computing  hardware. 

In  a  structural  sense,  ATD  lies  in  between  LES  and  DNS.  Large-scale  and  small-scale 
components  of  a  flow  field  are  treated  separately,  corresponding  respectively  to  the  resolved 
and  subgrid-scale  (SGS)  parts  of  LES,  but  (unlike  LES)  no  formal  averaging  or  filtering  is 
employ^  in  constructing  the  equations  of  ATD.  In  particular,  the  governing  equations  are 
split  into  large-  and  small-scale  parts  by  an  additive  decomposition.  There  arc  cross-terms 
connecting  these  equations  but,  in  contrast  to  LES,  there  are  equations  for  each  separate 
factor  in  these  cross-terms,  so  there  is  no  closure  problem.  An  additional  advantage  of  ATD 
is  that  the  algorithm  lends  itself  to  a  high  level  of  parallelization.  Thus,  on  computers  with 
many  parallel  processors,  the  clock  time  required  to  perform  a  simulation  can  be  greatly 
reduced.  Section  2.2  below  describes  a  parallel  ATD  algorithm. 

In  the  process  of  splitting  the  governing  equations  into  equations  for  the  large-  and 
small-scale  parts  of  a  flow,  ATD  requires  the  introduction  of  decomposition  parameters. 
These  may  be  selected  arbitrarily,  but  one  selection  at  least  has  a  clear  physical 
interpretation.  This  is  the  particular  decomposition  identified  by  Hylin  in  McDonough,  et 
al.(2],  in  which  the  split  equations  for  the  velocity  time-derivatives  are  in  the  form  of 
transport  equations  for  the  large-  and  small-scale  momentum.  For  incompressible  flow  these 
equations  are: 
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where  5  is  a  decomposition  divergence,  which  we  set  equal  to  zero.  These  equations  form 
the  basis  for  the  work  below,  except  for  the  1-D  work  on  chemical  kinetics  and  turbulence, 
and  the  work  described  in  the  appendix. 

1.2  Summary  of  the  Work  Performed  from  June  through  September,  1990. 

During  the  first  four  months  of  the  current  grant,  we  completed  the  initial  development  of  a 
numerical  algorithm  for  solving  the  small-scale  Navier-Stokes  equations  in  a  2-D  cell.  As  a 
test  problem,  we  imposed  periodic  boundary  conditions  on  the  cell  and  looked  at  the 
evolution  of  a  vortical  perturbation  under  the  influence  of  a  uniform  large-scale  shear.  The 
preliminary  results  demonstrated  that  the  algorithm  was  accurate  to  at  least  second-order  in 
time,  and  that  the  flow  seemed  to  be  related  to  other  flows  describable  in  terms  of  “inner 
coordinates.” 

Another  part  of  our  work  in  the  first  four  months  was  concerned  with  analyses  of  the 
2-D  implementation  of  ATD  in  generalized  coordinates.  The  goal  of  this  work  is  to  gain  an 
understanding  of  what  is  required  in  order  to  construct  reasonably  portable  ATD  modules 
that  can  be  added  to  existing  CFD  codes  to  enhance  their  turbulence  simulation  capabilities. 
We  chose  to  employ  a  Fourier-Galerkin  representation  of  the  small-scale  equations  and 
looked  at  the  issue  of  including  metric  information  in  the  Galerkin  inner  products.  We  also 
began  the  development  of  a  code  for  the  large-scale  Navier-Stokes  equations,  employing 
generalized  coordinates  in  two  dimensions,  llie  code  is  based  on  a  projection  method  like 
that  employed  by  Kim  and  Moin  [3],  extended  to  generalized  coordinates  and  using 
intermediate  boundary  conditions  that  result  in  second-order  temporal  accuracy.  The 
individual  split  steps  are  quasi-linearized  in  5-form  and  solved  via  Douglas-Gunn  time- 
splitting  [4J. 

Lastly,  we  continued  to  examine  the  mathematical  relations  between  ATD  and  related 
methods,  including  the  development  of  fast  solvers  for  elliptic  problems  using  domain 
decomposition  techniques. 


2.  TRANSITION  TO  TURBULENCE  IN  CHANNEL  FLOW 

2.1  Implementation  of  the  Small-Scale  Equations 

2.1.1  Summary  of  Work  to  Date.  In  our  most  recent  report  [5],  we  mentioned  our 
implementation  of  the  2-D  small-scale  equations.  This  uses  a  Chebyshev  tau  method  for  the 
spatial  discretization  of  the  equations,  the  Crank-Nicolson  method  for  the  temporal 
discretization  of  ±e  linear  velocity  terms,  and  a  second-order  predictor-corrector  method  for 
the  temporal  discretization  ot  the  nonlinear  and  pressure  gradient  terms.  The  resulting 
equations  are  in  Helmholtz  or  Poisson  form,  and  we  solve  them  by  means  of  a  variation  on 
the  Haidvogel-Zang  algorithm  [6,7]. 

We  first  implemented  the  small-scale  equations  using  periodic  boundary  conditions  on 
the  velocity  but,  as  discussed  below,  these  later  proved  unsuita*'!:  for  our  test  picblc.m. 
Furthermore,  a  full  implementation  of  the  ATD  algorithm  requires  a  combination  of  Dirichlet 
and  Neumann  conditions  on  the  boundaries  of  the  small-sc^e  cells.  Accordingly,  we  have 
modified  our  computer  code  for  the  small-scale  equations  to  handle  an  arbitrary  combination 
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of  Dirichlet  and  Neumann  conditions.  The  code  supports  up  to  33  Chebyshev  modes  in  each 
direction,  and  can  use  for  the  velocity  any  initial  condition  which  meets  the  mode  restriction 
and  is  zero  on  the  boundaries.  With  minor  modifications,  the  small-scale  code  is  ready  to  be 
used  as  a  subroutine  in  an  implementation  of  the  complete  ATD  algorithm  for  2-D  flow  in 
rectangular  geometries. 

2.12  Description  of  the  Test  Problem.  As  a  test  problem,  we  have  examined  the  decay  of  a 
vortical  perturbation  in  a  single  small-scale  cell,  subject  to  a  uniform  large-scale  shear.  We 
initially  imposed  periodic  boundary  conditions  [5],  but  imposition  of  the  proper  pressure 
boundary  condition  —  that  identifi^  by  Gresho  and  Sani  [8]  —  showed  that  the  small-scale 
pressure  was  not  periodic,  and  hence  that  periodic  boundary  conditions  were  not  appropriate. 
The  symmetry  inherent  in  the  definition  of  the  test  problem,  together  with  the  results 
obtained  using  periodic  boundary  conditions,  then  led  us  to  examine  rotationally  symmetric 
boundary  conditions  of  order  two.  Interestingly,  imposing  these  proved  to  be  identical  to 
imposing  homogeneous  Dirichlet  conditions  on  dl  four  boundaries,  i.e.,  to  using  a 
completely  closed  boundary,  which  is  also  inappropriate  for  this  flow. 

This  identity,  however,  led  to  a  test  of  a  projection  technique  which  will  be  used  in  an 
implementation  of  the  full  ATD  algorithm  to  ensure  that  the  small-scale  solution  is  globally 
continuous.  Briefly,  this  technique  uses  Lagrange  multipliers  to  project  onto  a  divergence- 
free  subspace  a  trial  solution  for  the  velocity  fluxes  through  the  boundaries.  We  us^  this 
technique  in  conjunction  with  the  condition  on  rotational  symmetry,  and  compared  the 
resulting  solution  to  that  obtained  by  imposing  homogeneous  Dirichlet  conditions  directly. 
The  two  solutions  were  identical  almost  to  machine  precision,  completely  validating  the 
projection  technique. 

The  appropriate  boundary  conditions  for  this  flow  turn  out  to  be  those  corresponding  to 
an  isolated  vortex  —  one  far  from  any  boundarier  to  the  large-scale  flow  —  or  a  “semi- 
isolated”  vortex  —  one  bounded  by  a  wall  on  one  side  only.  For  the  isolated  vonex,  open 
(i.e.,  flow-through)  boundary  conditions  are  used  on  all  four  sides  of  the  cell,  while  for  the 
semi-isolated  vortex,  open  conditions  are  used  on  three  sides  and  a  homogeneous  Dirichlet 
condition  is  used  on  the  side  representing  the  wall.  We  are  currently  modeling  open 
boundaries  by  imposing  an  homogeneous  Neumann  condition  on  the  velocities,  and  on  all 
boundaries  have  used  a  homogeneous  Neumann  condition  for  the  pressure.  This  latter 
condition,  especially,  is  not  quite  correct  [8],  but  a  recent  paper  by  Gresho  suggests  that  the 
errors  may  be  confined  close  to  the  bound^es  [9]. 

The  small-scale  test  problem  is  defined  in  a  square,  x  e  [-1,1],  y  e  [-1,1],  and  the 
initial  small-scale  perturbation  is  given  by  the  stream  function: 

V  =  C(l-4x^-)-6x'‘-4x®-i-x*)(  l-4yS6y'’-4y'^-i-y*)  (2-1) 

where  C  is  a  constant  expressing  the  strength  of  the  perturbation.  This  stream  function 
yields  a  vortex  which  has  zero  velocity  on  the  boundary  of  the  defining  square.  In  fact,  the 
velocity  is  on  the  boundary.  The  perturbation  pressure  is  determined  by  using  the 
penurbation  velocity  in  the  Poisson  equation  for  the  small-scale  pressure: 

Ap*  =  V- 

On  SOHO  boundaries  we  used  for  the  perturbation  pressure  the  Neumann  condition  identified 
bv  Gresho  and  Sani  [?];  with  the  given  stream  fiir,ctic.i,  this  becomes  an  hom^'jenenus 
Neumann  condition.  On  open  boundaries  we  set  the  perturbation  pressure  to  zero,  which 
implicitly  sets  the  initial  acceleration  on  those  boundmes.  Setting  the  initial  pressure  on 
open  boundaries  to  zero  also  ensures  that  the  pressure  is  C°  on  the  boundary.  Thus  defined. 
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both  the  velocity  perturbation  and  the  pressure  perturbation  have  compact  support. 

The  form  of  the  initial  vortex  is  given  by  equation  2-1.  If,  additionally,  the  initial 
strength  is  held  constant,  then  the  test  problem  is  governed  by  a  single  parameter,  which  may 
be  expressed  either  as  a  Reynolds  number  or  as  a  non-dimensional  vortex  diameter.  The  two 
are  related  according  to 


xV  ^  =  2VrJ 

'  V 


(2-3) 


where  X  is  the  diameter  of  the  initial  vortex,  du/dy  is  the  shear  strain-rate  for  the  large-scale 
flow,  and  Re  is  the  half-diameter  Reynolds  number. 


2.1  J  Results  to  Date  for  the  Test  Problem.  We  have  solved  the  problem  described  above  for 
Reynolds  numbers  ranging  from  1  to  3000.  For  the  cases  we  have  simulated,  figures  2-1  and 
2-2  show  the  decay  or  growth  of  the  Chebyshev-weighted  \}  norm  of  the  small-scale  velocity 
field: 


(2-4) 
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Figure  2-1.  Decay  of  the  velocity  norm  of  an  isolated  vortex. 


In  most  of  the  cases  shown,  there  is  an  initial  transient  —  due  to  our  particular  initial 
and  boundary  conditions  —  followed  by  a  long  period  of  decay.  The  rate  of  decay  becomes 
slower  and  more  unsteady  as  the  Reynolds  number  is  increased,  and  at  high  enough  Reynolds 
numbers  the  vortex  grows  in  strength.  Similar  behavior  is  shown  in  figure  2-3,  which 
focuses  on  the  growtn  or  decay  of  the  (0,1)  component  of  the  small-scale  x-velocity.  This 
component  is  constant  in  the  x-direction,  and  linear  in  the  y-direction,  and  for  low  Reynolds 
numbers  is  a  mode  which  decays  slowly.  For  Reynolds  numbers  below  10,  the  decay  of  the 
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(0.1)  mode  is  evident  in  the  figure.  For  Reynolds  numbers  trom  10  to  300,  the  initial 
transient  lasts  longer  than  the  time  frame  shown  in  the  figure,  but  we  believe  that  the 
corresponding  curves  do  evenmally  decay.  For  Reynolds  numbers  of  1000  and  3000, 
however,  an  initial  decay  is  followed  by  the  appearance  of  an  oscillation  which  grows  in 
amplitude.  Thus,  figure  2-3  suggests  a  transition  Reynolds  number  somewhere  in  the  range 
of  300  to  1000.  We  should  mention  here  that  the  time-step  size  for  these  simulations  was 
10^,  so  the  period  of  the  oscillations  is  many  times  the  time-step  size  and  thus  the 
oscillations  are  not  due  to  a  numerical  instability. 


Figure  2-2.  Decay  of  the  velocity  norm  of  a  semi-isolated  vonex. 


Figure  2-3.  Decay  of  the  (0,1)  mode  of  the  x-velocity  for  an  isolated  vonex. 
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Recall  that  the  Reynolds  number  in  this  problem  is  related  by  equation  2-3  to  a  non- 
dimensional  diameter  X*.  In  figure  2-4,  the  initial  and  late-stage  decay  rates  of  the  velocity 
norm  for  an  isolated  vortex  are  plotted  as  functions  of  X*.  A  late-stage  curve  for  the  semi- 
isolated  vortex  was  not  plotted  because  some  of  the  curves  in  figure  2-2  do  not  seem  to  have 
settled  down  to  an  asymptotic  decay  rate  by  the  end  of  the  simulation,  but  the  early-stage 
curves  for  the  semi-isolated  and  isolated  vortices  are  identical.  The  values  of  X"  which  are 
of  special  interest  are  those  where  the  curves  cross  the  zero  line  and  the  -1  line.  Vortical 
disturbances  smaller  than  the  former  value  tend  to  be  damped  out,  while  larger  disturbances 
are  amplified  and  may  begin  to  affect  the  large-scale  flow.  The  two  curves  in  figure  2-4 
cross  the  -1  line  at  X*  =  3  and  X*  =  8,  and  the  late  curve  seems  to  cross  the  zero  line  for  some 
X*  between  10  and  70.  Vortical  disturbances  with  X*  smaller  than  3  to  8  are  damped  so 
rapidly  that  they  have  a  negligible  effect  on  the  overall  flow;  they  decay  ir  less  time  than  it 
takes  the  fluid  nearby  to  travel  a  distance  comparable  to  the  vonex  diameter.  Vortical 
disturbances  with  X*  larger  than  a  transition  value  between  10  and  70  begin  to  grow  in 
strength.  The  results  in  figure  2-3  suggest  that  the  transition  value  lies  between 
X*  =  35  (Re=300)  and  X*  =  63  (Re=1000).  These  values  are  consistent  with  our  previous 
results,  even  though  the  boundary  conditions  are  different. 


X 

1  10  100 


Figure  2-4.  Growth  rate  of  the  velocity  norm  vs.  X* 


As  described  in  our  most  recent  report  [5],  X*  is  expressed  in  “inner  coordinates.”  The 
same  values  of  X*  declared  above  to  be  significant  are  also  significant  in  the  environment  for 
which  inner  coordinates  v/erc  originally  defined:  the  inner  layers  of  a  wall-bounded  turbulent 
flow.  In  such  a  situation,  X*  <5  corresponds  to  the  laminar  sublayer,  5  <  X*  <  70 
corresponds  to  the  transition  region,  and  X*  >  70  corresponds  to  the  fully  turbulent  outer 
flow.  The  close  correspondence  between  these  values  and  the  results  we  have  obtained  show 
that  the  latter  are  consistent  with  and  provide  a  possible  explanation  for  a  well-known  body 
of  experimental  data. 
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Figure  2-5.  Evolution  of  stream  function  contours  for  an  isolated  vortex. 


The  structural  evolution  of  the  vortex  is  affected  both  by  the  large-scale  shear  and  by 
the  influence  of  the  Reynolds  number  (or  the  non-dimensional  vortex  iameter).  Figure  2-5 
shows  the  evolution  of  the  stream  function  contours  for  an  isolated  vortex  at  a  Reynolds 
number  of  300.  The  shear  flow  moves  to  the  right  at  the  top  of  the  cell  and  to  the  left  at  the 
bottom,  and  the  vortex  is  rotating  counter-cloclwise.  As  might  be  expected,  the  large-scale 


shear  causes  a  tilting  and  elongation  of  the  vortex,  but  it  also  tends  to  cause  an  intensification 
of  the  vorticity,  an  effect  which  becomes  apparent  at  high  Reynolds  numbers.  This  is 
illustrated  in  figure  2-6,  which  shows  the  vorticity  contours  for  three  cases; 
Re  =  100  (X*  =  20),  Re  =  300  (V  =  35),  and  Re  =  1000  (X'*'  =  63),  as  well  as  for  the  initial 
condition.  When  Rc  =  100,  the  vonicity  gradually  decays,  but  when  Re  =  300  or  above,  the 
vorticity  becomes  more  intense  with  time.  This  suggests  that  transition  for  diis  flow  may 
begin  in  the  neighborhood  of  Re  =  300. 


Figure  2-6.  Intensification  of  vorticity  at  high  Reynolds  numbers. 


The  structural  evolution  of  a  semi-isolated  vortex  is  qualitatively  similar  to  the 
evolution  of  an  isolated  vortex,  but  —  in  the  case  when  the  vortex  is  rotating  counter¬ 
clockwise  —  the  presence  of  the  wall  induces  negative  vorticity  and  positive  flow  between 
the  vortex  and  the  wall.  This  is  evident  in  figure  2-7,  where  the  wall  is  on  the  bottom  side  of 
each  square.  The  decay  of  a  semi-isolated  vortex  is  prolonged  by  the  presence  of  a  virtual 
vortex  on  the  other  side  of  the  wall  but,  as  shown  by  a  comparison  of  figiu'es  2-1  and  2-2,  the 
presence  of  the  wall  also  causes  qualitative  changes  in  the  decay  transient.  At  present,  the 
reasons  for  these  changes  are  not  known. 
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Str.  Fn,,  Re  «  300,  T  -  0  5,  Bottom  Wall  b.  Stream  Fr...  Re  »  300,  Time  -  1.5 


c.  Vorticity,  Re  -  300,  Time  •  0.5  d.  Vortioty,  Re  -  300,  Time»i  .5 


Figure  2-7.  Evolution  of  a  senii-isolated  vortex. 


While  figures  2-1  through  2-3  and  rigure  2-6  show  results  for  Re  =1000  and 
Re  =  3000,  and  while  the  early  results  for  these  cases  look  acceptable,  we  were  unable  to 
carry  the  integrations  for  these  cases  beyond  a  certain  point  without  the  computations 
blowing  up.  In  these  cases,  however,  we  have  observed  that  the  vortex  appeared  to  be 
stretched  beyond  the  confines  of  the  simulation  cell.  It  may  be  that  as  the  vortex  starts  to 
grow,  it  begins  to  exert  a  significant  influence  on  the  flow  beyond  the  cell  boundaries,  which 
is  rot  included  in  our  model.  This  should  not  be  a  pioblem  in  the  complete  ATD  algorithm 
(except  possibly  at  outflow  boundaries),  because  the  neighboring  cells  will  then  be  present 
and  flow  between  cells  will  be  accounted  for.  At  outflow  boundaries  the  boundary 
c'nditions  will  have  to  be  treated  carefully;  certainly  more  carefully  than  we  have  treated 
them  in  this  test  problem. 

2.2  Definition  of  a  Complete,  Parallelizabie  ATL  Algorithm 

We  have  continued  to  proceed  with  the  implementation  of  the  complete  ATD  procedure  for 
the  Navier-Stokes  equations,  applied  to  the  case  of  2-D  channel  flow,  and  have  identified  an 
algorithm  which  in  some  circumstances  can  '•e  completely  parallelized.  This  algcnthm 
incorporates  and  makes  extensive  use  of  the  previously  dercribed  algorithm  for  the  2-D 
small-scale  equations  [5],  In  the  text  below,  we  shall  first  present  the  algorithm  and  then 
discuss  '^ich  step  in  detail.  In  the  listing  of  the  algorithm,  each  level  of  indentation  indicates 
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either  a  loop  structure  or  parallelization  at  the  subroutine  level. 


Algorithm 

Predict  the  large-s<"ale  solution  at  the  next  large-scale  time  step. 

For  each  small-scale  cell: 

Project  the  large-scale  solution  onto  the  small-scale  basis. 

For  each  small-scale  time  step: 

For  each  small-scale  cell: 

Interpolate  the  large-scale  solution  for  the  current  time. 

Carry  out  the  preliminary  small-scale  solution. 

While  not  converged: 

For  each  cell  interface: 

Smooth  the  velocities  across  the  interface. 

.Project  the  smoothed  velocities  onto  a  divergence-free  subspace. 
For  each  cell: 

Resolve  the  small-scale  problem. 

For  each  small-scale  ^ell: 

Restrict  the  solution  onto  the  large-scale  bast'. 

Correct  the  large-scale  solution  at  the  new  large-scale  time. 


2.2.1  Predict  the  large-scale  solution  at  the  next  large-scale  time  step.  This  is  the  first  step 
in  the  ATD  solution  algorithm  defined  by  McDonough  [2,10].  Supposing  that  the  over-all 
algDiithm  is  to  be  p*  order  accurate,  and  that  n  large-scale  time  steps  have  been  completed. 
In  this  step  of  the  current  algorithm,  we  calculate  a  p-1*  order  prediction  of  the  large-scalw 
solution  at  the  n-i-1*  step.  We  will  probably  implement  this  step  serially,  using  the  same 
solution  routine  we  have  developed  for  the  small-scale  equations.  However,  as  long  as  the 
time-discretization  is  carried  out  such  that  the  resulting  semi-discrete  equations  are  in 
Poisson  or  Helmholtz  form,  the  many  techniques  whi*  h  have  been  developed  for  the  parallel 
solution  of  Poisson's  equation  may  be  used  to  construct  a  parallelized  procedure  for  this  step 
in  the  algorithm,  .\ssuming  that  we  use  the  serial  routine  that  we  have  on  hand,  if  there  are  N 
small-scale  cells  in  an  approximately  square  array,  and  one  large-scale  data  point  per  cell, 
this  step  will  require  approximately  O  (N^^)operations. 

2.22  ..Froject  the  large-scale  solution  onto  the  small  scale  basis.  After  transmission  of  the 
large-scale  solution  to  each  processor,  this  step  can  be  carried  out  in  parallel  on  as  many 
processors  as  there  are  small  scale  cells.  If  there  are  M  small-scale  data  points  in  each  cell 
(M'^)  modes  in  each  direction,  this  step  will  take  0(Nx(M’^N  -t- MN‘^ -i- M  log  M  ) ) 
operations. 

2.23  For  each  small-scale  time  step...  The  small-scale  time  step  is  ^pically  much  smaller 
than  the  large-scale  time  step.  Tnis  statement  in  the  algorithm  begins  a  loop  which  runs  until 
the  solution  for  the  small-scale  flow  has  evolved  from  the  n*  large-scale  timt  step  to  the 

n+l‘". 

2.2.4  Interpolate  the  large-scale  solution  for  the  current  time.  The  statement:  "For  each 
small-scale  cell:"  begins  a  section  of  code  which  is  parallelizable.  The  first  step  in  this 
section  is  to  interpolate  the  large-scale  solution  for  the  current  small-scale  time.  Recall  that 
the  large-scale  solution  has  been  estimated  to  p-l*  order  accuracy  at  the  n-»-l*  large-scale 
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time.  In  order  for  the  small-scale  solution  to  be  accurate,  and  because  the  small-scale  time 
step  is  much  smaller  than  the  large-scale  time  step,  interpolation  of  the  large-scale  solution 
between  the  values  at  the  n*  and  n-t-l'*'  times  is  required.  If  p  =  2,  a  linear  interpolation  will 
suffice,  in  which  case  this  step  requires  0(M)  operations  per  cell,  or  0(NM)  operations 
altogether. 

2.2  J  Carry  out  the  preliminary  small-scale  solution.  Once  the  large-scale  solution  has  been 
interpolated  to  the  current  small-scale  time,  a  preliminaiy  solution  of  the  small-scale 
equations  is  carried  out  in  each  cell.  This  may  be  done  in  parallel  if  there  is  more  than  one 
processor.  Using  our  small-scale  solution  routine,  this  step  requires  operations  per 

cell,  or  O  (NM^^)  operations  altogether.  The  results  are  only  preliminary  because  we  do  not 
yet  have  sufficient  information  to  accurately  update  the  velocities  on  the  cell  boundaries. 
Getting  this  information  requires  iteration. 

2.2.6  While  not  converged...  This  statement  begins  a  loop  in  which  the  velocities  on  the  cell 
boundaries  are  updated  and  the  solution  to  the  small-scale  equations  is  recalculated.  This  is 
repeated  until  the  solution  converges.  The  convergence  criterion  is  a  function  of  the  large- 
scale  and  small-scale  time  steps,  as  well  as  of  the  accuracy  required  on  the  large  scale. 

2.2.7  Smooth  the  velocities  across  the  interface.  In  any  solution  of  an  elliptic  problem  which 
is  carried  out  by  means  of  a  domain  decomposition,  it  is  necessary  to  smooth  the  solution 
across  the  boundaries  of  the  subdomains.  TTiis  may  be  accomplished  by  various  methods, 
several  of  which  we  have  already  examined  in  connection  with  ATD  [10,11,12].  Here  we 
propose  to  fit  a  fifth-order  polynomial  across  each  boundary,  matching  the  small-scale 
velocity  solutions  with  continuity  at  the  nearest  Chebyshev  collocation  point  on  either 
side.  Evaluating  the  polynomial  at  the  boundary  will  then  provide  an  updated  value  for  the 
boundary  velocity. 


Figure  2-8.  Smoothing  of  the  small-scale  velocity  profile. 


This  technique  is  shown  schematically  in  figure  2-8,  where  the  solid  lines  represent  the  un¬ 
smoothed  solution  and  the  dashed  lines  represent  the  polynomial.  Because  the  fitted 
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polynomial  blends  smoothly  wiiii  the  old  velocity  profiles,  the  new  boundary  value  for  the 
velocity  at  this  stage  satisfies  the  small-scale  momentum  equations  to  at  least  O  (h),  where  h 
is  the  distance  from  the  boundary  to  the  nearest  collocation  point  This  step  in  the  algorithm 
requires  O  ( M  log  M )  operations  per  cell,  or  0(  NM  log  M )  operations  adtogether.  It  may 
be  parallelized,  in  which  case  communication  is  required  between  the  processors  holding 
data  for  adjacent  cells. 

2.2.8  Project  the  smoothed  velocities  onto  a  divergence-free  subspace.  This  is  a  crucial  step, 
and  was  alluded  to  in  section  2.1.  The  smooth^  velocities  on  the  cell  boundaries  do  not 
necessarily  satisfy  the  requirement  for  mass  conservation.  In  this  step,  the  mass  flux  across 
each  cell  boundary  is  first  computed  using  the  new  boundary  velocities.  Then,  Lagrange 
multipliers  are  employed  to  adjust  the  fluxes  until  they  satisfy  the  zero-divergence  (mass 
conservation)  constraint,  while  changing  the  fluxes  as  little  as  possible.  The  constraints  on 
this  optimization  problem  do  not  change  —  they  are  determined  by  geometry  alone  —  but 
the  data  (the  cell  divergences)  are  different  each  time  the  problem  is  solved.  The  problem 
can  be  reduced  to  a  set  of  matrix/vector  multiplications,  requiring  O(N^)  operations,  and 
there  are  many  existing  algorithms  which  may  be  used  to  parallelize  this  step. 

2.2.9  For  each  cell,  resolve  the  small-scale  problem.  In  each  cell,  the  small-scale  equations 
are  solved  once  again,  using  the  updated  values  for  the  velocities  on  the  boundaries.  Note 
that  only  the  boundary  values  for  the  small-scale  velocities  are  updated.  The  boundary 
condition  for  the  pressure  is  always  determined  in  each  cell  from  the  momentum  equation 
applied  on  the  boundary  (c.f.  Gresho  and  Sani  [8]).  As  before,  this  step  requires  0(M^^) 
operations  per  cell,  or  O  (NM^'^)  operations  altogether. 

2.2.10  Restrict  the  solution  onto  the  largescale  basis.  Once  the  small-scale  solution  has 
been  advanced  to  the  n-f-l*  large-scale  time  step,  it  is  necessary  to  restrict  the  small-scale 
solution  onto  the  large-scale  basis,  so  that  the  large-scale  solution  may  be  corrected.  This 
requires  O  (NM)  operations  altogether,  and  the  restriction  may  be  carried  out  in  all  cells  in 
parallel. 

2.2.11  Correct  the  large-scale  solution  at  the  new  large-scale  time.  This  is  the  last  step  in 
the  ATD  solution  algorithm  identified  by  McDonough,  and  corrects  the  previous  large-scale 
solution  to  p*  order  at  the  new  time.  Once  again,  if  we  use  the  same  code  we  have  been 
using  to  solve  the  small-scale  equations,  this  step  will  require  approximately 
O  (r'I^^)operations. 

2.2.12  Summary  of  operation  counts.  If  we  take  K  time  steps  on  the  small  scale  for  each  time 
step  we  take  on  the  large  scale,  and  if  L  iterations  are  required  at  each  small-scale  step  for  the 
solution  to  converge,  then  one  step  through  the  above  algorithm  requires 

O  (M^^^  -H  MN^'^  -H  MN  log  M  KM^'^  -i-  KLM^'^  -K  KLN^  -t-  KLMN  log  M) 

operations.  Because  the  algorithm  is  completely  parallelizable,  many  of  these  operations  can 
be  carried  out  simultateously  if  more  than  one  processor  is  available.  Minimizing  execution 
clock-time,  however,  will  depend  on  minimizing  the  amount  of  communication  required 
between  processors,  as  well  as  the  operation  count  of  each  parallelized  step.  This  will  impel 
careful  consideration  of  how  the  data  are  distributed  to  the  various  processors. 
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3.  ATD  FOR  2-0  INCOMPRESSIBLE  FLOW  IN  GENERALIZED  COORDINATES 


Our  current  work  on  ATD  in  generalized  coordinates  is  divided  into  two  parts:  first, 
implementation  of  the  large-scale  equations  using  a  chaotic  map  to  simulate  the  effects  due  to 
the  small  scale;  and  second,  implementation  of  the  small-scale  equations  with  an  imposed 
large-scale  flow.  Not  until  both  implementations  have  been  checked  out  and  understood  will 
they  be  combined  to  form  a  complete  ATD  procedure. 

3.1  The  Large-Scale  Equations  in  Generalized  Coordinates 

The  code  that  we  are  developing  to  solve  the  large-scale  incompressible  Navier-Srokes 
equations  in  2-D  generalized  coordinates  is  based  on  a  projection  method.  In  such  a  method, 
the  Navier-Stokes  momentum  equations  are  solved  with  the  pressure-gradient  term  omitted, 
and  the  solution  is  then  projected  onto  a  divergence-free  velocity  field.  The  resulting 
solution  satisfies  the  continuity  equation  exactly,  while  the  momentum  equation  is  satisfied 
only  approximately,  but  the  error  in  the  momentum  equation  can  be  kept  small  enough  for 
the  method  to  be  of  practical  use. 

The  derivation  of  the  large-scale  equations  in  generalized  coordinates  is  complex  and 
is  here  omitted,  as  are  the  equations  themselves.  Details  are  in  Appendix  A.  As  is  generally 
the  case  in  ATD,  there  are  cross-terms  in  the  large-scale  equations  which  involve  the  small- 
scale  velocities.  Initially,  instead  of  obtaining  the  small-scale  velocities  by  solving  the  small- 
scale  equations,  we  will  model  them  by  means  of  a  chaotic  map  which  we  will  develop,  with 
the  expectation  that  the  large-scale  equations  derived  via  ATD  and  used  in  conjunction  with 
such  a  map  may  be  able  to  simulate  more  general  turbulent  flows  than  is  possible  with  LES 
and  the  current  models  for  subgrid-scale  effects. 

So  far  we  have  derived  the  generalized-coordinate  versions  of  the  large-scale 
momentum  equations  without  the  pressure  gradient  terms.  These  equations  have  been 
implemented  in  VS  FORTRAN  on  the  CBM  3090-600J  supercomputer  at  UK,  and  constitute 
the  first  step  in  the  projection  algorithm.  Our  next  tasks  in  this  area  are  to  implement  the 
projection  step,  and  to  develop  the  chaotic  map  to  model  the  small-scale  velocities. 

3.2  The  Small-Scale  Equations  in  Generalized  Coordinates 

We  have  decided  to  use  a  Fourier-Galerkin  procedure  for  the  solution  of  the  small-scale 
equations  in  generalized  coordinates,  in  part  because  we  believe  that  it  will  be  applicable  for 
a  wide  range  of  problems,  and  in  part  because  it  is  easier  to  implement  in  the  context  of 
generalized  coor^nates  than  is  the  Chebyshev-tau  approach  used  in  the  channel  flow 
problem.  In  particular,  it  is  jxjssible  to  arrange  things  such  that  the  basis  functions  satisfy  the 
small-scale  continuity  equation  (1-4)  mode  by  mode.  This  requires  the  introduction  of  the 
contravarient  velocities  U  and  V: 

U  =  ^^u-i-^v  y  =  T]u  +  r{\,  (3-1) 

where  ^ ,  etc.  are  elements  of  the  Jacobian  matrix  of  the  coordinate  transformation. 
Then  the  smil-scale  continuity  equation  in  generalized  coordinates  becomes 


which  has  a  structure  very  similar  to  that  in  the  Canesian  coordinates.  Here  J  =  . 

To  construct  the  local  approximation,  we  consider  a  cell  centered  at  (4, ,  q  j  and 
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assume  that  the  solution  can  be  represented  by  truncated  Fourier  series,  defined  as  follows. 


(3-3a) 

(3-3b) 

(3-3c) 

where 

=  kjc,  Oj  =  l;t. 


•  .  Tl-Tl 

U  (^,Ti,t)  =  J  2^a(t)yCosa^— J^sino,— J— 

ki 

V' 

V  (^,Ti,t)  =  J  ^by(t)sin(x^— ^cosOj — j— 

kj 

P  (^,ri.t)  =  2jCu(0sina^  -— since, — j- — 


and  h  is  the  finite  difference  grid  spacing  for  the  large-scale  calculation.  A  major  advantage 
of  in  choosing  the  basis  functions  as  defined  by  equation  3-3  is  that  when  we  construct 
Galerldn  inner  products  we  can  eliminate  the  Fourier  coefficient  for  pressure  in  a  way  similar 
to  that  suggested  by  Orszag  and  Kruskal  [13]  (See  also  Canuto,  et  al,  [7]). 


(3-4) 


where  f^  and  g^  are  the  Galerkin  inner  products  of  the  nonlinear  terms.  (See  Appendix  A). 

Using  the  basis  functions  defined  in  equation  3-3,  we  have  derived  the  Fourier- 
Galerkin  representation  of  the  small-scale  equations  in  generalized  coordinates.  The 
derivation  of  the  small-scale  momentum  equations  in  generalized  coordinates  is  complex  and 
is  here  omitted,  as  are  the  equations  themselves  —  details  are  in  Appendix  A  —  but  the 
derivation  gives  a  set  of  ODE's  for  the  coefficients  a^i  and  by.  We  solve  this  system  of  ODE's 
via  a  second-order  explicit  Runge-Kutta  method  (Heun's  method).  We  will  be  using  the  code 
based  on  these  equations  to  investigate  the  convergence  of  Fourier-Galerkin  approximations 
to  chaotic  solutions  of  the  Navier-Stokes  equations,  and  to  investigate  ^e  effects  of 
bifurcation  parameters  such  as  the  Reynolds  number  and  the  local  gradients  of  the  large-scale 
velocity  on  the  nattire  of  solutions  to  the  small-scale  Navier-Stokes  equations. 


4.  CHEMICAL  KINETICS/TURBULENCE  INTERACTION 

It  is  extremely  difficult  to  simulate  the  details  of  interactions  between  turbulence  and  other 
physical  phenomena.  Reynolds-averaged  approaches,  and  LES,  employ  modeling  on 
precisely  the  length  and  time  scales  where  interactions  are  likely  to  be  important,  and  DNS  is 
too  expensive  to  be  able  to  produce  the  desired  results  for  complex  engineering  problems. 
During  the  past  year  we  have  begun  to  study  the  use  of  the  small-scale  equations  of  ATD  as 
local  parametrized  (with  large-scale  properties)  DNS.  We  have  completed  an  initial 
simulation  corresponding  to  the  tip  of  an  H2-O2  diffusion  flame  at  low  Mach  number  (- 
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0.05),  and  moderate  Re  (-  6000).  Details  oi  these  studies  are  reported  by  McDonough  and 
Saito  [14]. 

Here,  we  begin  by  presenting  the  governing  equations  being  used,  in  dimensionless 
form,  the  1-D,  (nearly)  compressible,  viscous  Navier-Stokes  and  species  equations: 

p,  +  (pu),  =  0  (4- la) 

(pu)^  +  (pu"\  -  =  -P,  (4-lb) 

T.  +  T„  =  -^(PY,)^  (PYJ  e-'"*  (4.10 

(pY,)+(puY)  -^(PYJ  =-C^-f(pY,/(pY,)e-'"’^,  i=  1.  2.  (4-ld,c) 

t  X  i  Cq  X 

In  these  equations 


0,  ,0, 
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where  Pex  and  Peo  are,  respectively,  thermal  and  mass  diffusion  Peclet  numbers; 


As 
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(4-3) 


and 


^L(P°T  J  B 

uVjWj 


■-(P°T  J  B 
uV^ 


(4-4) 


is  a  reference  temperature  (different  from  'P’),  and  B  is  the  pre-exponential  factor.  T  is 
the  dimensional,  unshifted  temperature,  and  t  is  the  corresponding  dimensionless  quantity. 

We  apply  the  additive  turbulent  decomposition  in  the  manner  described  by 
McDonough  and  Saito  [14],  and  then  using  the  Fourier  representation 


K 

p*(x,t)  =  ^a^(t)sina^x*  ,  (4-5a) 

k=l 

K 

m*(x,t)  =  ^b|^(t)cosa^x*  ,  (4-5b) 

k=l 
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(4-5c) 


K 

T*(x,t)  =  ^c^(t)smOj^x*  , 

k=l 

K 

z*(x,t)  =  ^d^(t)sinO|^x* , 

k=l 

K 

z*(x,t)  =  ^e^(t)sinaj^x*  , 

k=l 

where  =  kic/h,  and  x*  s  x-(x -h/2),  with  xe  rx-h/2,  x.+h/2l ,  we  construct 
Galerkin  procedure  which  leads  to  the  small-scale  ODE’s 


i  =a  b 

^  kh  k 


ijJ=l 


(4-5d) 

(4-5e) 
a  local 

(4-6a) 

(4-6b) 

(4-6c) 

(4-6d) 

(4-6e) 

(4-7) 
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K. 

+  ^  rz,dc,+2z,c.c.+td.e.')d,c  F... 

V  2  I  j  I  I  J  1  J/  I  m  yli 

+  ^  d  d  e.c  .c  G  ,  ’e  , 

m  n  I  J  1  ijlrnnk 


fVT 

Vik=J  ^sina^: 

f'T 

cosa^ 

'iik=  cosa^: 

•'«.4 


x*sinaj^x*sinay,x*  dx 


x*sinOj^x*coscru,x*  dx 


x*cosap,x*cosay,x*  dx 


f'-T 


;^x*coscTjhX*sinOy,x*  dx 


ri 

s  sinaihX*sir 

«'v| 

fVT 

F>iimk  =  sinai,x*si 


sincyjhX*sinaii,x*sinOu,x*  dx 


sinaj^x*sina„,x*sina^x*sinay,x*  dx 


=  sina^x*sinajhX*sinay.x*si 


sina,^x*sina^x*sinay,x*  dx 


(4-8a) 


(4-8b) 


(4-80 


(4-8d) 


(4-8e) 


(4-8f) 


(4-8g) 


In  addition,  we  have  the  trivial  system  for  k  =  1; 
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Figure  4-1  Small-Scale  Fluctuations  vs.  Time:  a)  velocity,  b)  temperature, 
c)  H2  concentration,  d)  O2  concentration 


The  Galerkin  ODE’s,  Eqs.  (4-6),  are  nonlinear  and  for  large  K  can  be  expected  to  be 
stiff.  Thus,  one  would  typically  expect  to  employ  Gear’s  method  as  the  solution  technique. 
But  because  we  are  anticipating  chaotic  solutions  (and  hence,  sensitivity  to  initial 
conditions),  and  because  Gear’s  method  is  a  multi-step  method  (with  potential  for  parasitic 
error),  we  have  chosen  to  employ  a  modification  of  the  A-stable,  single-step  trapezoidal 
scheme.  The  modification  consists  of  treating  the  nonlinear  terms  explicitly  (via  Heun’s 
method),  while  maintaining  implicit  treatment  of  the  linear  terms,  ^e  main  source  of 
stiffness  in  this  case.  We  of  course  do  lose  unconditional  stability,  but  we  also  gain  a 
significant  reduction  in  the  arithmetic  that  would  be  required  to  solve  the  linear  systems  that 
would  have  arisen  in  the  application  of  Newton’s  method  to  the  nonlinear  terms  had  they 
been  handled  implicitly.  This,  in  turn,  reduces  the  round  off  error,  which  is  highly  desirable 
when  integrating  dynamical  systems  possessing  chaotic  solutions,  again  because  of 
sensitivity  to  initial  conditions. 

Seven  modes  were  retained  in  the  small-scale  representations  of  all  dependent 
variables;  i.e.,  K  =  7.  The  dimensionless  parameter  values  were  Re  =  6000,  Pex  =  ICKX)  and 
Peo  =  3000.  Notice  that  we  have  not  employed  a  Le  =  1  assumption.  The  physical  size  of 
our  small-scale  region  in  the  neighborhood  of  the  flame  tip  was  ().5  mm.  The  dimensionless 
time  step  size  employed  in  our  numerical  integrator  was  At  =  5.0x10^,  corresponding  to  a 
physical  time  step  of  0.666...  |isec.  The  simulations  were  run  for  100  time  steps. 

Results  for  velocity,  temperature  and  H2  and  O2  concentration  fluctuations  are 
displayed  in  Fig.  4-1.  These  were  obtained  by  inserting  the  solutions  to  Eqs.  (4-6)  at  each 
time  t  into  Eqs.  (4-5)  and  summing  only  from  ic  =  2  to  k  =  K.  This  constmcticn  starts  with  k 
=  2  because  k  =  1  is  actually  a  Fourier  mode  corresponding  to  pait  ot  the  large-scale  solution. 

There  are  a  number  of  interesting  observations  to  be  made  regarding  Fig.  4-1.  To 
begin,  the  temporal  behavior  on  the  small  scale  corresponds  to  what  should  be  termed 
transient  chaos.  Careful  examination  of  the  figures  shows  that  oscillatory  behavior  first 
begins  in  the  small-scale  velocity  and  temperature  almost  simultaneously.  But  this  occurs 
only  after  a  relatively  long  period  of  increase  in  negative  amplitude  of  the  H2  and  O2  small- 
scale  concentrations,  which  is  consistent  with  depletion  of  reactants  to  form  products.  From 
Eq.  (4-6b)  we  see  that  as  the  negative  small-scde  concentration  amplitudes  become  large, 
they  lead  to  positive  forcing  of  the  momentum  equation,  provided  the  large-scale  local 
temperature  gradient  is  positive.  We  have  estimated  this  to  have  a  value  of  825  K/cm  from 
data  of  Saito  et  al.  [15].  Once  the  oscillations  in  velocity  and  temperamre  begin,  they  feed 
back  into  the  species  equations  through  the  advective  terms  (for  velocity),  and  the  species 
production  terms  (for  temperature.).  It  can  also  be  seen  that  negative  concentration 
fluctuations  result  in  generally  positive  temperature  fluctuations,  as  would  be  expected  on 
physical  grounds,  since  the  negative  concentration  fluctuations  imply  product  formation,  and 
thus  heat  release.  In  general,  these  results  appear  to  be  generally  consistent  with  expected 
physics. 
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5.  MATHEMATICAL  STUDIES  RELATING  TO  ATD 


Mathematical  studies  of  techniques  related  to  ATD  were  conducted  at  UCLA  and  focused  on 
two  issues.  One  focus  was  on  a  comparison  of  a  technique  for  coupling  the  large  scale  and 
small  scales  in  the  2D  incompressible  Navier-Stokes  equations,  to  obtain  efficient  long  time 
integrations.  The  other  focus  was  on  the  development  of  fast  solvers  for  elliptic  problems  in 
two  dimensions,  using  domain  decomposition  techniques.  Such  solvers  can  be  used  to  solve 
the  linear  systems  arising  from  semi-implicit  time  stepping. 

The  original  version  of  the  ATD  method  was  based  on  a  decomposition  of  the  solution 
into  two  scales,  a  large  scale  and  a  small  scale  component.  Each  component  of  the  solution 
was  computed  by  an  iterative  procedure.  The  studies  at  UCLA  incorporated  the  recently 
developed  non-linear  Galerkin  method  proposed  by  Temam  [16],  Marion  and  Temam  [17] 
and  Titi  [18],  which  was  also  based  on  a  decomposition  of  the  solution  into  two  sedes. 
However,  the  two  scales  were  coupled  through  a  stationary  map  which  defined  the  small 
scale  component  in  terms  of  the  large  scale  component,  in  such  a  way  that  the  resulting 
solution  was  on  a  non-linear  manifold  approximating  the  global  attractor  of  the  problem.  The 
domain  decomposition  studies  focused  on  a  comparison  of  this  technique  for  coupling  the 
two  scales  with  the  standard  pseudo- spectral  method  in  which  the  two  scales  were  coupled 
cxticily. 

Tests  conducted  for  the  2D  incompressible  Navier-Stokes  equations  indicated  that  the 
non-linear  Galerkin  method  can  cost  less  than  a  corresponding  pseudo-spectral  method,  due 
to  the  possibility  of  choosing  larger  time  steps.  The  coupling  used  in  [16]  between  the  two 
scales  resulted  in  good  approximations,  when  the  number  of  modes  used  to  represent  the 
large  scale  component  was  sufficiently  large.  However,  the  method  developed  some 
instabilties  when  the  number  of  large  scale  males  was  not  sufficiently  large.  This  restriction 
corresponded  to  the  number  of  modes  required  to  resolve  the  flow,  as  predicted  by  the  theory 
of  Henshaw,  Kreiss  and  Reyna  [19].  More  tests  are  needed  to  study  the  possibilty  of  more 
efficient  hybrid  algorithms  which  adapts  the  number  of  modes  used  as  the  flow  evolves  into 
large  scale  structures. 

The  studies  at  UCLA  on  the  development  of  fast  solvers  for  elliptic  problems,  focused 
on  using  domain  decomposition  techniques  with  non-overlapping  subdomains.  Efficient 
versions  of  algorithms  were  developed  having  rates  of  convergence  independent  of  mesh 
parameters,  and  also  of  coefficient  variations,  for  symmetric  positive  definite  elliptic 
problems  in  2  dimensions.  In  some  cases,  the  algorithms  had  a  complexity  of  0(N  log  N), 
where  N  denotes  the  number  of  unknowns. 

Details  of  these  domain  decomposition  studies  are  presented  in  Appendix  B. 


6.  SUMMARY  OF  WORK  TO  DATE 

The  work  performed  under  the  two  grants  —  AFOSR  90-0271  and  89-0281  —  has  covered  a 
range  of  subjects  related  to  Additive  Turbulent  Decomposition.  Work  with  Burgers’ 
equation  resulted  in  the  development  of  a  non-iterative  algorithm  for  coupling  the  results  of 
large-scale  and  small-scale  time  integrations,  and  the  rigorous  proof  of  its  consistency  and 
accuracy.  These  results  were  subsequentiy  extended  to  the  Navier-Stokes  equations.  A 
straightforward  physical  argument  was  found  which  determines  a  unique  set  of  vdues  for  the 
decomposition  tensors  required  when  ATD  is  applied  to  the  Navier-Stokes  equations.  A 
study  of  the  effects  of  Reynolds  averaging  on  the  solution  of  Burgers’  equation  demonstrated 
that  the  ensemble  average  of  the  solution  to  the  unaveraged  equations  equals  the  solution  to 
ihe  ensemble  averaged  equations  if  and  only  if  properly  averaged  Reynolds  stresses  are  used 
in  the  latter,  and  that  the  time  average  of  the  solution  to  the  unaveraged  equations  is  never 
exactly  equal  to  the  solution  of  the  time  averaged  equations,  even  when  exact  Reynolds 
stresses  arc  employed. 

The  2-D  small-scale  equations  have  been  implemented  using  a  Chebyshev-tau  method. 
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and  applied  used  to  examine  the  decay  or  growth  of  an  isolated  vortex.  This  problem  is 
related  to  the  origin  of  turbulence  in  -hear  flows.  The  results  suggest  tha*  this  problem  is  also 
related  to  the  structure  of  wall-bounded  flows,  and  tend  to  confirm  the  la^ge  body  of 
experimental  data  associated  with  such  flows.  The  code  developed  to  examine  this  pi  blem 
will  be  at  the  heart  of  a  larger  code  implementing  a  .highly  parallel  ATD  algorithm  we  have 
defined  for  the  Navier-Stokes  e'luations.  Work  has  begun  on  the  application  of  ATD  to  the 
incompressible  Navier-Stokes  equations  i-i  2-D  generalized  coordinates.  The  large-scale  and 
small-scale  equations,  expressed  in  ge.'.’ralized  coordinates,  have  been  derived,  and  their 
implementation  is  in  progress. 

A  number  of  mathematical  studies  have  been  carried  out  regarding  ATD  and  domain 
decomposition.  Two  particular  domain  decomposition  methods  were  sullied  in  connection 
with  ATD  applied  to  Burgers’  equation.  The  results  showed  that  the  Schwarz  alternating 
method  was  more  robust  than  the  Gauss-Seidel-Newion  (GSN)  metnod,  and  required  fewer 
iterations,  but  the  arithmetic  per  iteration  was  much  lower  for  the  GSN  method  than  for  the 
Schwarz  method.  An  approach  related  to  ATD,  the  nonlinear  Galerkin  method  of  Temam 
[16],  Marion  and  Ternam  [17]  and  Titi  [18],  was  compared  to  the  pseudo-specual  method  of 
Henshaw,  Kreiss  and  Reyna  [19],  and  proved  to  be  somewhat  more  efficient.  And  studies 
have  been  made  of  efficient  and  easily  parallelizable  methods  for  solving  elliptic  problems. 
These  are  especially  applicable  to  the  solution  of  the  small-scale  equations  on  the  global 
domain. 

Future  work  is  expected  to  include  the  implementation  of  the  complete  ATD  algorithm 
for  the  Navier-Stokes  equations  and  its  application  to  the  problem  of  transition  in  channel 
flow,  solution  of  the  large-scale  equations  using  a  chaotic  map  to  model  the  subgrid-scale 
effects,  a  parametric  examination  of  the  bifurcations  in  the  solutions  to  the  small-scale 
equations,  and  the  development  of  a  complete  ATD  code  for  generalized  coordinates. 
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APPENDIX  A:  Derivation  of  the  ATD  Equations  in  Generalized  Coordinates 


APPENDIX  A:  DERIVATION  OF  THE  ATD  EQUATIONS  IN  GENERALIZED 
COORDINATES 

A.l  Derivation  of  the  Large-Scale  Momentum  Equations  in  Generalized  Coordinates 

With  the  pressure-gradient  terms  omitted,  the  Navier-Stokes  momentum  equations  in  2-D 
cartesian  coordinates  are: 


3u  3  ,  2,  3  1 


ay 


Re 


Define  a  general  coordinate  transformation  T  such  that 


(A-1) 

(A-2) 

(A-3) 


(  Xj  Xj  )<-T^(  ^  ) 

We  can  apply  this  transformation  to  equations  A-1  and  A-2  and,  using  the  chain  rule  for 
derivatives,  obtain  the  corresponding  equations  in  generalized  coordinates: 


and 
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Using  the  particular  decomposition  described  by  Hylin  in  McDonough,  et  al.  [A-1], 
large-scale  momentum  equations  corresponding  to  equations  A-1  and  A-2  are 


+  -^(u^  +  -^(vu)  +  -^(0*0)  -I-  ■^(v*u)  =  ^Au 
ot  ox  ay  ox  ay  Re 

-I-  -^(uv)  +  +  -^(u*v)  -t-  •^(v*v)  =  -^Av 

ot  dx  ay  dx  ay  Re 
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where  u*  and  v*  denote  the  small-scale  velocities.  In  generalized  coordinates,  these 
equations  become: 
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The  underscored  terms  are  those  involving  the  small-scale  velocities  u*  and  v*,  and  do  not 
appear  in  equations  A-4  and  A-5.  We  will  be  developing  a  chaotic  map  to  model  u*  and  v*, 
with  the  expectation  that  the  large-scale  equations  deriv^  via  ATD  and  used  in  conjunction 
with  such  a  map  may  be  able  to  simulate  more  general  turbulent  flows  than  is  possible  with 
LES  and  the  current  models  for  subgrid-scale  effects. 

In  discretizing  equations  A- 8  and  A-9,  we  first  linearize  them  by  means  of  the  5-form 
quasilinearization,  then  rearrange  the  equations  to  make  the  diffusive  and  convective  terms 
clearer.  For  the  spatial  discretization  we  use  the  standard  staggered  grid,  with  centered 
differencing  for  the  diffusive  terms  and  first-order  upwinding  for  the  convective  terms. 

A.2  Derivation  of  the  Small-Scale  Momentum  Equations  in  Generalized  Coordinates 

As  derived  by  Hylin  in  McDonough,  et  al.  [A-1],  the  small-scale  momentum  equation  is 

=  -^Au*  -  V-(uu*)  -  V-(u*u*)  -  Vp*  (A- 10) 

dt  Re 

while  in  generalized  coordinates,  tlie  momenmm  equations  are  in  the  form 

Q-t-F.-t-G  =-;^(R„-f-S,  +T  ),  (A-11) 


as  given  by  Fletcher  [A-2].  The  contravariant  velocities  U  and  V  are  defined  to  be 

U  =  ^^u-(-^v  V  =  q^u TiyV  ,  (A-12) 


where  ,  fly ,  etc.  are  elements  of  the  Jacobian  matrix  of  the  coordinate  transformation.  The 
momentum  equations  contain  both  contravariant  velocities  and  primitive  velocities  and  in 
terms  of  these  velocities  the  terms  in  the  small-scale  version  of  equation  A-11  become 
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where  U  and  V  are  the  large-scale  contravariant  velocities. 

To  construct  the  local  approximation,  we  consider  a  cell  centered  at  ,  ti.),  and 
assume  the  solution  can  be  represented  by  truncated  Fourier  series,  defined  as  follows. 
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where 


=  kTi,  a,  =  In, 

and  h  is  the  finite  difference  grid  spacing  for  the  large-scale  calculation.  We  now  substitute 
Eq.(A-14)  into  (A- 13)  and  construct  the  Galerkin  inner  products.  The  resulting  system  of 
ordinary  ifferential  equations  for  the  time-dependent  Fourier  coefficients  is 
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where 


C,  =  J Jsin^-i^  =  JJcos^^  dtidri , 

and  fy  and  are  the  Galerkin  inner  products  of  the  nonlinear  terms: 
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1  Introduction 

This  project  concerns  a  study  of  some  techniques  for  efficiently  integrating  the  incom¬ 
pressible  Navier-Stokes  equations  over  large  time  intervals.  For  such  problems,  with 
time  independent  forcing,  turbulent  solutions  may  converge  in  function  space  to  a 
finite  dimensional  set  referred  to  as  the  global  attractor,  see  [4].  In  this  case,  the  long 
term  behaviour  of  the  solutions  can  essentially  be  described  by  a  finite  dimensional 
sytem  of  ODE’s  which  captures  the  large  scale  properties  of  the  global  attractor. 
Thus  far,  however,  computing  accurate  solutions  of  the  incompressible  Navier-Stokes 
equations  by  direct  simulation  is  still  very  expensive,  especially  in  three  dimensions. 
Some  alternate  techniques  have  been  proposed  which  reduces  this  cost,  such  as  Large 
Eddy  Simulation  techniques  (LES),  see  [5]  [6].  Recently,  the  .Additive  Turbulent  De¬ 
composition  method  (ATD)  was  proposed  by  McDonough  and  Bywater  [9]  [10]  as  an 
alternative  direct  simulation  technique  without  the  use  of  statistical  parameters.  In 
this  report,  we  breifly  describe  our  studies  on  the  ATD  method  and  related  meth¬ 
ods  for  numerically  integrating  the  incompressible  Navier-Stokes  equations  over  large 
intervals  of  time. 

We  compare  two  approaches.  In  one  case,  we  follow  a  method  related  to  the  basic 
idea  in  the  original  version  of  the  ATD  method  proposed  in  [9,  10],  which  was,  to 
decompose  the  solution  as  a  sum  of  a  large  scale  and  small  scale  component  and 
to  approximately  determine  each  component  by  an  iterative  technique  based  on  the 
original  equations.  Our  studies  led  us  to  incorporate  the  coupling  used  in  the  recently 
developed  non-linear  Galerkin  method  of  Temam  [13],  Marion  and  Temam  [8]  and 
Titi  [14].  The  non-linear  Galerkin  method  is  also  based  on  a  splitting  of  the  solution 
into  two  components,  large  scale  and  small  scale.  However,  unlike  the  coupling  used  in 
the  original  .ATD  version,  the  coupling  between  the  large  scale  and  small  scale  terms 
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is  determined  by  a  stationary  mapping  which  defines  the  small  scale  component  in 
terms  of  the  large  scale  component.  We  use  this  form  of  coupling  between  the  two 
scales. 

In  another  approach,  we  use  a  pseudo-spectral  method  of  Henshaw,  Kreiss  and 
Reyna  [7],  in  which  the  number  of  modes  used  is  chosen  to  provide  a  convergent 
numerical  approximation.  This  is  based  on  a  theory  for  the  minimum  scale  of  the 
flow,  presented  in  [7].  The  minimum  scale  can  also  be  used  to  adapt  the  number  of 
modes  required  to  resolve  the  flow,  as  the  flow  evolves  into  coherent  structures,  therby 
leading  to  improved  efficiency.  In  Section  5,  we  compare  the  non-linear  Galerkin 
method  with  the  pseudo-spectraJ  method.  Our  studies  are  restricted  to  the  vorticity- 
stream  function  formulation  of  the  2D  incompressible  Navier-Stokes  equations  on  a 
periodic  domain. 

The  rest  of  the  report  is  outlined  as  follows.  We  describe  the  vorticity-stream  func¬ 
tion  formulation  of  the  incompressible  Navier-Stokes  equations  on  a  periodic  domain 
in  2  dimensions.  We  then  describe  its  discretisation  by  a  pseudo-spectral  method 
in  space  and  a  4th  order  predictor  corrector  in  time.  Following  that  we  describe 
the  non-linear  Galerkin  method.  We  present  some  numerical  comparisons  of  the  two 
methods  considered  and  some  of  their  properties  are  discussed.  Finally,  we  report  on 
some  related  work  we  have  done  on  developing  fast  solvers  for  elliptic  problems  using 
domain  decomposition  techniques.  Such  algorithms  can  be  applied  to  solve  the  linear 
systems  obtained  when  semi-implicit  time  stepping  is  used  for  non-periodic  boundary 
conditions. 

2  Navier-Stokes  equations  in  two  dimensions 

The  velocity-pressure  formulation  of  the  2  dimensional  incompressible  Navier-Stokes 
equations  on  the  square  [0,2^]  x  [0,2t]  with  periodic  boundary  conditions  reads: 

U(  +  uu^  +  vUy  +  Vp  =  t'Au -t- f  in  [0,27r]2 

V  •  u  =0  in  [0, 27r]2  (1) 

u(t  =  0, .)  =  Uq  for  t  =  0 

where  u=  (u,  v)  denotes  the  velocity,  p  denotes  the  pressure,  f  denotes  a  time  inde¬ 
pendent  forcing,  and  ly  denotes  the  kinematic  viscosity.  The  vorticity  tu  =  V  x  u 
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satisfies: 

UJt  +  UCJ^  +VUJy  =  l/U  +  f 

Uj{t  =  0,x,y)  =  uJoix,y)  ^  > 

where  (u,u)  =  u  =  {tpy,  —4'x)  's  the  curl  of  the  stream  function  ly  which  satisfies; 

—Alp  =  ijj. 

Thus,  in  the  2  dimensional  case,  u,  v  and  rp  are  determined  by  the  vorticity  and  we 
will  consider  the  solution  of  scalar  evolution  equation  (2)  for  u). 

3  A  pseudo-spectral  method 

We  now  describe  a  pseudo-spectral  method  used  in  Henshaw,  Kreiss  and  Reyna  [7], 
and  Browning  and  Kreiss  [2],  which  leads  to  convergent  numerical  approximations 
for  the  vorticity  equation  (2).  In  the  pseudo  spectral  method,  the  solution  a;  is 
approximated  by  a  truncated  Fourier  series  Uy,  whose  coefficients  (which  depend  on 
time)  evolves  as  a  system  of  ODE’s: 

/V/2-1 

ui{t,x,y)  ^  A{t,  k^, 

fcl,*2  =  -,V/2+l 

where,  for  each  k=  Fourier  coefficient  d;(t,  ^1,^:2)  satisfies: 

k)  +  k)  -I-  vuPy{t.  k)  -f  i/\k\'^u!{t,  k)  =  /(k). 

The  ODE  system  can  be  written: 

iDj  -f  ui\w  =  G(th), 

where  G{w)  =  —uWj.  —  vWy  +  f  represents  the  convection  term  and  the  forcing,  and 
A  is  a  diagonal  matrix  with  entries;  =  |kp. 

In  our  tests,  the  ODE  system  for  the  coefficients  u;(k)  were  integrated  by  the  use 
of  a  4th  order  Predictor-Corrector  method,  eis  in  [7],  in  which  the  diffusion  term  was 
integrated  exactly.  For  time  step  At  the  prediction  step  reads: 

lip  =  -f-  ^  (23e-‘'^^‘G„  -  , 

and  the  correction  step; 

-f  ^  (9Gp  +  ■ 
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In  our  tests,  for  mesh  patrameter  h  =  2/N,  the  time  step  At  was  chosen  to  satisfy 
the  following  CFL  condition: 

<  (|uL  +  kloo)  ^ 

h 

where  CFL^i^  chosen  to  be  O.i  and  CFL^^j.  was  chosen  to  be  1.2.  So  long  as 
the  choice  of  At  kept  the  CFL  quantity  v/ithin  the  bounds,  the  step  size  was  not 
altered.  But  when  it  fell  outside  the  range,  the  new  time  step  was  chosen  to  be: 


0.5h 

“  (IxU  +  H*)’ 

The  cost  of  evaluating  the  derivitive  u,\  is  roughly  proportional  to  the  cost  of 
computing  the  convolution  term.  This  can  be  evaluated  approximately  using  5  two 
dimensional  FFT’s  by  mapping  cJ...  u  and  v  to  real  space,  forming  uu^  +  vu:^  in 
real  space  and  mapping  the  result  back  to  Fourier  space  (pseudo-spectral  method). 

The  theory  of  Henshaw,  Kreiss  and  Reyna  [7]  gives  an  explicit  bound  for  the 
maximum  number  of  m'  des  needed  to  resolve  the  flow: 


where  v  is  the  kinematic  viscosity  and  is  the  maximum  of  the  velocity  gradients 

for  the  flow,  for  all  time,  if  that  is  known  a  priori.  Over  large  time  intervals,  as  the  flow 
evolves  into  large  scale  structures,  it  would  be  possib'e  to  adapt  for  improved 

efficiency,  by  using  the  above  estimate.  We  intend  on  testing  the  efficiency  of  adaptive 
gridding. 


4  Non-linear  Galerkin  method 

As  mentioned  earlier,  our  studies  of  the  coupling  between  the  two  scales  in  the  original 
version  of  the  ATD  method,  led  us  to  consider  an  alternate  form  of  coupling  used 
in  the  non-linear  Galerkin  method  of  Marion  and  Temam  [8]  and  Titi  [14].  We  have 
used  thij  coupling  in  the  tests  reported  here.  We  now  briefly  describe  the  nonlinear 
Galerkin  method  for  an  evolution  equation  in  a  function  space  H: 

-f  oAw  -I-  B{w)  =  f  . 

u;(0)  =Wo 


4 


where  y4  denotes  a  self-adjoint  positive  operator  on  H  (such  as  -A),  u  denotes  a  small 
parameter,  B{w)  is  a  nonlinecir  map  (such  as  uw,, -)- utWy),  and  /  iS  a  time  independent 
forcing.  The  function  space  is  decomiposed  H  =  where  represents 

the  space  spanned  by  ''ome  large  scale  eigenfunctions  of  A  and  represents  its 
orthogonal  complement.  The  solution  w  can  be  decomposed  w  =  p  +  q,  where  p  € 
and  q  G  .  Equation  (3)  is  then  equivalent  to  the  coupled  system: 

Pt^  uAp+ P^Bip  +  q)  =  , 

qt  +  i/Aq  +  PsB{p  +  q)  =  Psf 

where  P[^  denotes  orthogonal  projection  onto  and  Ps  =  I  —  Pi-  In  the  original 
version  of  the  ATD  method,  system  (4)  was  solved  by  means  of  an  iterative  procedure. 
However,  in  the  non-linear  Galerkin  method,  undor  the  assumption  that  ||(^||  IIpII 
and  <C  II^A^-r  Ps(B(p)  —  /)||,  the  2nd  equation  in  (-t)  is  appro.ximated  by: 

i^Aq  +  PsB{p)  ~  Psf, 

which  definv,s  the  small  scale  component  q  in  terms  of  the  large  scale  component  p, 
by  a  function  ^(p)  defined  by 

q  =  ^{p)  =  {i\A)~^Ps{f  -  B{p)).  (5) 

The  resulting  new  system  for  the  large  scale  component  p  reads. 

Pt  +  uAp  +  PiB{p  +  ^{p))  =  P^f  (6) 

Equation  (6j  comprises  the  nonline.ar  Galerkin  method  for  p,  and  is  different  from  the 
standard  Galerkin  equations  for  p:  t  includes  the  nonlinear  coupling  between  the  two 
scales.  For  each  p,  its  corresponding  small  scale  component  ^{p)  can  be  determined 
using  (5),  thereby  leading  to  a  solution  in  Hi  -f  Hs- 

In  our  application,  the  space  H  C  Z/^([0,2:r]  x  [0. ‘Jt|)  corresponds  to  the  Fourier 
modes  of  wavenumbers  less  than  iV,  Hi  corresponds  to  Fourier  modes  of  wavenumbers 
less  than  .V/2  and  Hs  corresponds  to  the  remaining  Fourier  modes  in  H .  In  the  tests, 
we  discretised  evolution  equation  (6)  for  the  large  scale  component  v,  by  an  application 
of  the  pseudo-spectral  method.  We  note  that  the  method  described  here  can  also  be 
e.xtended  to  the  case  where  finite  difference  or  finite  element  discretisations  are  used, 
as  hcis  been  reported  in  the  1991  ICI.AM  conference. 
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5  Numerical  Tests 


We  present  here  tlie  csnits  of  a  sample  test  comparing  the  performance  of  the  non¬ 
linear  Gah  rkui  method  with  the  standard  pseudo-spectral  method,  as  described  in 
Sections  3  and  4.  In  particular,  we  present  contour  plots  of  the  vorticity,  plots  of  the 
energy  “^p',  trum,  the  spectral  decay  in  the  vorticity,  and  plots  of  the  norms  of  the  laxge 
scale  and  mall  scale  components  of  the  solution,  for  the  two  methods  considered. 

In  the  ‘est,  we  chose  me  forcing  /  to  be  zero.  The  initial  data  ujq  for  the  vorticity 
was  chosen  so  :hat  the  Fourier  coefficients  ^>f  the  corresponding  stream  function  tl’o 
satisfied; 

,,  t,.|=/c.{lk|-^l!+(|k|/6f|}‘'h  |k|<20.5  1 
'  [0,  |k|  >  20.5  J 

as  in  Browining  am’  Kreiss  [2]  The  phases  of  the  Fourier  coefficients  were  chosen 
randomly  from  [0,27r].  .A.s  described  in  [2],  if  we  chose  Ci  so  that  =  1,  then 

we  can  appro.ximately  evaluate  |Z)u|,^.  For  this  particular  choice  of  C-,  and  for  ^ 
kinematic  viscosity  u  =  5.0  x  IQ—’,  we  obtain  an  estimate  for  the  nu;  iber  of  mo  c; 
required  to  resolve  the  flow  to  be: 


-v; 


34. 


In  the  test,  we  therefore  chose  64  modes  in  each  direction,  for  the  pseudo-spectral 
method,  with: 

|yti|<3l  |fc2l<3l 

For  the  non-linear  Galerkin  method,  we  chose  32  modes  in  e^-ch  direction  for  the  large 
scale  component  (extended  by  its  small  scale  component  so  that  there  are  still  64  x  64 
modes  in  their  sum),  i.e., 

^.VGAf  =  P  +  <7* 


where 


and 


P=  E  E 

|ki  |<15  16<1*:2|<31 

16<ifci|<31  |t2li;15  16<l*i|<31  16<l*2l<31 
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The  total  number  of  modes  for  p  -f  q  is  therefore  64  x  64. 

For  the  above  choice  of  initial  data  ljq,  the  energy  E{1)  is  proportional  to  l~^  for 
large  /,  where  E[l)  represents  the  energy  corresponding  to  wavenumbers  of  magnitude 
/: 

^  X]  (1^1^  +  1^1^)  =  T  for /,  =  [/- 0.5,/ +  0.5).  (7) 

“  ik|e/,  "  |k|6/, 

Similarly,  \  quantity  W (/)  representing  the  average  magnitude  of  the  Fourier  coeffi¬ 
cients  of  of  wavenumber  /  is  defined: 

W{1)  =  (  Y:  kk)|]  /  f  E  l]  •  (8) 

Vlk|€/,  /  \|k|e/, 

These  quantities  are  computed  in  the  numerical  tests. 

In  Figure  1,  we  present  the  contour  plots  of  the  the  vorticity,  obtained  by  both  the 
pseudo-spectral  method,  and  the  non-linear  Galerkin  method,  at  time  T  ~  5.0.  In 
Figure  2,  the  contour  plots  at  time  T  ^  40.0  is  presented.  The  times  are  not  exactly 
40.0,  since  variable  time  steps  were  chosen  by  using  the  CFL  condition.  In  Figure  3  we 
:)mpare  the  energy  spectrum  of  the  solution  at  time  T  «  5.0  and  T  w  40.0  obtained 
by  both  methods.  Figure  4  is  a  similar  plot  of  the  quantity  4F(/)  defined  in  (8),  at 
times  T  v  5.0  and  T  «  40.0.  In  Figure  5  we  plot  the  norms  of  the  components  llp||, 
llqll  and  ||q(||  as  functions  of  time,  to  see  whether  the  approximation  that  ||q||  <C  ||p|| 
and  llqjj  <t.  j|p||  is  accurate. 

A  discussion  on  the  numerical  results.  The  results  in  Figures  1  through  5 
indicate  that  the  non-linear  Galerkin  method  based  on  N/2  modes  (plus  the  remaining 
modes  defining  the  small  scale  component  q)  gives  solutions  which  are  similar  to  those 
obtained  by  the  pseudo-spectral  method  with  N  modes.  Since  only  half  the  number 
of  modes  were  used  for  p  in  equation  (6),  the  non-linear  Galerkin,  can  use  times  steps 
At  which  are  twice  as  large  as  that  required  by  CFL  limitations  for  the  corresponding 
pseudo-spectral  method  with  N  modes.  However,  the  evaluation  of  ^(p)  as  well  as 
j5(p+<I>(p))  requires  some  2D  FFT’s  involving  N  modes  in  each  direction,  rather  than 
N/2  modes  in  each  direction.  Overall,  there  is  a  reduction  in  computational  cost  by 
a  fixed  fraction,  due  to  the  presence  of  the  larger  time  step  allowed  in  the  non-linear 
Galerkin  method. 

rtlong  another  line,  it  was  observed  that  if  the  number  of  modes  N/2  used  to 
represent  the  large  scale  solution  p  was  not  sufficiently  large,  then  the  function  ^*(p) 


Figure  1;  Contour  plots  of  the  vorticity  at  T  ft:  5.0 


Result  for  pseudo-spectral  at  T=5.04 
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Result  for  non-linear  Galerkin  at  T=5.08 
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Figure  2:  Contour  plots  of  the  vorticity  at  I  ^  40.0 


Contour  interval  is  0.1,  Nl=n2=64,  nu=5.0E-4 
Result  for  non-linear  Galerkin  at  7=40.07 


Contour  interval  is  0.1,  N=64,  N/2=32,  nu=5.0E-4 


Figure  5:  Norms  of  !|p|i,  and  lj'7t||. 


xlO'^  Pseudo-spectral  method.  Norms  of  "p",  ”q"  and  "q_t "  ,vs.  time 


llpll  marked  "o";  llqli  marked  llq_tll  marked 
xlO-^  Non-linear  Galerkin  method.  Norms  of  "p”,  "q"  and  "q_t”  .vs.  time 
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beccime  large  and  an  instability  developed.  We  attribute  this  to  the  following.  Since, 
by  the  definition  of  $(p): 

4(k  )  =  $(p)(k)  =  (//IkP)”*  {/(k)  —  5(p)(k)  j ,  for  indices  |k|  >  cN 

it  follows  that  if  the  divisor  i/|kp  is  small,  then  the  corresponding  Fourier  modes  are 
magnified,  and  an  instability  can  result.  This  can  be  fixed  by  requiring  the  number 
of  modes  N/2  to  be  sufficiently  large,  so  that: 

i/jkp  >  1,  for  |k|  >  cN. 

This  requirement  leads  indirectly  to  a  similar  condition  appearing  in  the  minimum 
scale  derived  by  Henshaw,  Kreiss  and  Reyna  [7]: 

N  > 


for  some  positive  constant  C. 

We  intend  to  perform  more  tests  on  these  methods,  including  the  case  where  the 
number  of  modes  used  is  varied  as  the  flow  evolves  into  large  scale  structures,  and 
also  the  use  of  hybrid  algorithms.  These  may  lead  to  accurate  integrations  over  large 
time  intervals  with  even  more  reduced  computational  cost. 

6  Elliptic  solvers 

Finally,  we  briefly  summarise  results  of  related  studies  on  efficient  and  easily  par- 
allelisable  methods  for  solving  elliptic  problems  [3],  which  can  be  applied  to  those 
linear  systems  which  occur  in  time  stepping  the  Navier-Stokes  equations.  Our  studies 
have  been  based  on  two  well  known  domain  decomposition  algorithms  [1]  [12].  They 
are  based  on  a  partition  of  the  domain  into  many  non-overlapping  subregions,  and  a 
preconditioner  is  developed  for  the  resulting  reduced  interface  problem.  The  precondi¬ 
tioner  essentially  corresponds  to  a  block  Jacobi  method,  where  the  blocks  correspond 
to  couplings  between  unknowns  on  certain  subregions  of  the  interface.  The  subregions 
of  the  interface  are  chosen  to  be  the  edges  separating  the  subdomains,  a  coarse  grid 
consisting  of  the  vertices  of  the  subdomains,  and  cross-shaped  regions  centered  about 
the  vertices.  Our  studies  have  focused  on  replacing  the  exact  blocks  of  the  reduced 
interface  matrix,  by  spectrally  equivalent  preconditioners,  as  this  significantly  reduces 


13 


•i  %  •' 


the  over  head  cost.  Numerical  experiments  confirm  theoretical  results  indicating  that 
the  resulting  preconditioner  has  an  optimal  rate  of  convergence  [3].  With  the  use  of 
fast  Poisson  solvers  on  the  subdomains,  the  resulting  algorithm  can  be  implemented 
in  0(A^log(iV))  operations  where  N  is  the  number  of  unknowns.  Such  algorithms  can 
also  be  applied  to  indefinite  linear  systems,  such  as  those  obtcuned  by  discretisations 
of  the  velocity  -  pressure  formulation  of  the  Navier-Stokes  equations  [11]. 

References 

[1]  J.  H.  Bramble,  J.  E.  Pasciak  and  A.  H.  Schatz,  The  construction  of  precondi¬ 
tioners  for  elliptic  problems  by  substructuring,  /,  Math.  Comp.,  47(175):  103-134, 
1986. 

[2]  G.  L.  Browning  and  H.-O.  Kreiss,  Comparison  of  numerical  methods  for  the 
calculation  of  two-dimensional  turbulence,  Mathematics  of  Computation,  Vol. 
52,  No.  186,  pages  369-388,  April  1989. 

[3]  T.  F.  Chan  and  T.  P.  Mathew,  An  application  of  the  probing  technique  to  the 
vertex  space  method  in  domain  decomposition.  Tech.  Report  90-22,  UCLA,  De¬ 
partment  of  Mathematics,  Oct.  1990.  Submitted  to  Proceedings  of  Domain  De¬ 
composition  Conference  held  in  Moscow. 

[4]  P.  Constantin,  C.  Foias  and  R.  Temam,  Attractors  representing  turbulent  flows, 
in  Mem.  Amer.  Math.  Soc.  53,  No.  314,  1985. 

[5]  J.  H.  Ferziger,  in  Computational  methods  for  turbulent,  transonic  and  viscous 
flows,  J.  E.  Essers  (de.).  Hemisphere  Pub.  Washington,  1983. 

[6]  J.  H.  Ferziger,  in  Theoretical  approaches  to  turbulence,  Dwoyer  et  al  (eds.). 
Springer- Verlag,  New  York,  1985. 

[7]  W.  D.  Henshaw,  H.-O.  Kreiss  and  L.  G.  Reyna,  On  the  smallest  scale  for  the 
incompressible  Navier-Stokes  equations,  CAM  Report  88-10,  UCLA,  Depeirtment 
of  Mathematics,  March  1988. 

[8]  M.  Marion  and  R.  Temam,  Nonlinear  Galerkin  Methods,  SIAM  J.  Numer.  Anal., 
Vol.  26,  No.  5,  pp.  1139-1157,  October  1989. 


14 


[9]  J.  M.  McDonough  and  R.  J.  By  water,  Turbulent  solutions  from  an  unaveraged 
additive  decomposition  of  Burgers’  equation,  Forum  on  Turbulent  Flows  -  1989, 
Bower  and  Morris  (eds.),  ASME  FED  Vol.  76,  New  York,  1989. 

[10]  J.  M.  McDonough  and  R.  J.  Bywater,  Large  scale-effects  on  local  small-scale 
chaotic  solutions  to  Burgers’  equation,  AIAA  Journal,  Vol.  24,  No.  12,  pp.  1924- 
1930.  December  1986. 

[11]  J.  Pasciak,  Two  Domain  Decomposition  Techniques  for  Stokes  Problems,  In  T. 
Chan,  R.  Glowinski,  J.  Periaux  and  0.  Widlimd,  editors,  Second  International 
Symposium  on  Domain  Decomposition  Methods  for  Partial  Differential  Equa¬ 
tions,  Philadelphia,  PA,  1988,  SIAM. 

[12]  B.  F.  Smith,  An  optimal  domain  decomposition  preconditioner  for  the  finite  ele¬ 
ment  solution  of  linear  elasticity  problems.  Technical  Repot  482,  Courant  Insti¬ 
tute,  1989.  To  appear  in  SIAM  J.  Sci.  Stat.  Comput. 

[13]  R.  Temam,  The  Nonlinear  Galerkin  Method  in  Computational  Fluid  Dynamics, 
To  appear. 

[14]  E.  S.  Titi  On  approximate  inertial  manifolds  to  the  N avier-Stokes  equations,  in 
Journal  of  mathematical  analysis  and  applications,  Vol.  149,  No.  2,  July  1,  1990, 
Academic  Press,  New  York. 


15 


